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Density shock waves in confined microswimmers 
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Motile and driven particles confined in microfiuidic channels exhibit interesting emergent be¬ 
havior from propagating density bands to density shock waves. A deeper nnderstanding of the 
physical mechanisms responsible for these emergent structures is relevant to a number of physical 
and biomedical applications. Here, we study the formation of density shock waves in the context of 
an idealized model of microswimmers confined in a narrow channel and subject to a uniform external 
flow. Interestingly, these density shock waves exhibit a transition from ‘subsonic’ with compression 
at the back to ‘supersonic’ with compression at the front of the population as the intensity of the 
external flow increases. This behavior is the result of a non-trivial interplay between hydrodynamic 
interactions and geometric confinement, and is confirmed by a novel quasilinear wave model that 
properly captures the dependence of the shock formation on the external flow. These findings can 
be used to guide the development of novel mechanisms for controlling the emergent density distri¬ 
bution and average population speed, with potentially profound implications on various processes 
in industry and biotechnology such as the transport and sorting of cells in flow channels. 

PACS numbers: 47.63.Gd, 87.18.Hf, 05.65.-|-b 


The transport of micro-particles in flow channels occurs 
in a variety of industrial and natural processes. Examples 
include filtration [1] , colloid deposition [2] , droplet-based 
microfluidics [3] , blood cells in microflows [1] , sperm cells 
in the fallopian tube [5], and pathogens in the urinary 
tract [5]. The ability to control the density distribution 
and group velocity of micro-particles in flow would there¬ 
fore have numerous applications in physics and biology. 
This letter analyzes, via discrete particle simulations and 
macroscopic models, the emergent patterns in popula¬ 
tions of motile particles driven by an external flow in a 
rectangular micro-channel. 

The dynamics of single and many-swimmer systems 
are typically studied in unconfined settings [THg. The 
effects of confinement to narrow flow channels are less 
well understood; see, e.g., [Mil] and references therein. 
Recent experiments on driven droplets 0 (Mm and 
self-propelled colloids im I19j in quasi two-dimensional 
(Hele-Shaw type) channels show the emergence of trav¬ 
eling density waves, including density shocks at the wave 
front. Similar observations are reported in simulations 
of confined self-propelled particles |20|, albeit with den¬ 
sity shocks forming at the back of the wave. The pri¬ 
mary mechanism responsible for the emergence of these 
density waves is attributed to hydrodynamic interactions 
(HI) among the confined particles [HI [20] • Confinement 
in Hele-Shaw type geometries induces a distinct hydrody¬ 
namic signature; the particle far-held disturbance is that 
of a source dipole irrespective of the details of the particle 
transport mechanism, driven or self-propelled, pusher or 
puller 0II11IM5I1I12]. Additional conhnement to a nar¬ 
row channel does not alter the nature of hydrodynamic 
interactions but imposes impenetrability conditions on 
the lateral boundaries of the channel. These effects are 
properly accounted for in the model used in |20j . and. 
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FIG. 1. (a) Schematic of microswimmers in a Hele-Shaw rect¬ 
angular channel, (b) Dipolar far-field flows induced by the 
self-propelled motion and the background flow, (c) Sidewall 
confinement is accounted for using an infinite image system. 


thus, do not explain the discrepancy in the location of the 
shock formation between the simulations |20j and exper¬ 
iments [ISKIH]. Here, we show that the background flow 
and geometric confinement coupled with the HI among 
the swimmers and their motility properties lead to a rich 
variety of emergent behaviors, consistent with experi¬ 
mental evidence but unobserved in simulations before. 
Most notably, we show the emergence of density shocks 
that transition from ‘subsonic’ to ‘supersonic’ as the in¬ 
tensity of the background flow increases. We conclude 
by demonstrating that this richness can be exploited to 
control the density distribution and collective speed of 
the microswimmers. 

We briefly review the model of N microswimmers in 
Hele-Shaw confinement presented in [23] [24] based on the 
work in Eiiiia . Namely, one has 

Zn = f7e“*“" -I- pLw{zn) + Vn, Un = Re [z^wie'""] (1) 

Here, denotes the swimmer’s position in the complex 
plane z = x iy {i = V—1) and a„ its orientation, 
n = The overline (•) and the operator Re[-] 
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FIG. 2. Density shock waves in confined microswimmers. Left: (a) for V < c, compression shock wave forms at the back; (b) 
for V = c, shock is suppressed; (c) for V > c, compression wave forms at the front. The ‘speed of sound’ of the medium ahead is 
c = U + fiV. Middle: the density p is the local area fraction computed as in [3]. Simulations results of the quasi Id continuum 
model are superimposed in red. Right: the probability distribution of swimming speeds in the discrete simulations are shown in 
solid lines, and the speed of sound of the medium ahead is indicated in dashed lines. Parameter values are: W = 60, $a = 0.3, 
i/ = 1, snapshots taken at t — 1000. 


denote the complex conjugate and the real part, respec¬ 
tively. Each swimmer has a self-propelled speed U and is 
advected by the local flow velocity w{z). All variables are 
non-dimensional using the characteristic length and time 
scales dictated by the size and speed of the individual 
swimmer, p, and v are translational and rotational mo¬ 
bility coefficients whose values depend on the geometric 
properties of the swimmer and its flagellar activity: p is 
in the range 0 < /r < 1 while v is either positive or nega¬ 
tive depending on the head-tail asymmetry [22] or flagel¬ 
lar activity [23]. We use a collision avoidance mechanism 
Vn based on the repulsive part of the Leonard-Jones po¬ 
tential. These near-held interactions decay rapidly out¬ 
side a small excluded area centered around thus en¬ 
suring that the order of the far-held HI is preserved. 

To close the model in Q, we need to evaluate the 
velocity held w{z) created by N conhned swimmers in 
uniform background how, (Fig. &)■ The self-propelled 
motion of a swimmer zj induces a dipolar velocity held 
Ws{z — Zj) = a^U/{z — ZjY which points in the swim¬ 
ming direction, whereas the background how induces an¬ 
other dipolar held Wb{z — zj) = —a^{l — p)Vl{z — Zj)'^ 
which points in the direction opposite to the background 
how, irrespective to the orientation of the swimmer (see 
Fig. & and Ref. [H]). Here, a and U are the effective 
radius and self-propelled speed of the swimmer, normal¬ 
ized to 1, and V is the velocity of the background how. 
This model is similar to the one used in [2nj . However, 
the fact that Wb is proportional to (1 —/r) overlooked 
in [20], thus failing to correctly capture the dipolar held 
induced by the background how. The velocity held w{z) 
created by N swimmers in uniform background how can 
be written as a superposition of Ws and Wb- 

The impenetrability conditions at z = ±iIF/2 due 
to the sidewall conhnement are accounted for using 
the method of images. Basically, to each swimmer at 
Zi = Xi + ij/i correspond two periodic lattices, of pe- 


V<c 


> 0: shock at front 


tr, < 0: shock at back 



FIG. 3. Id model does not capture correctly the shock wave 
formation. It results in compression shock at the front for 
V < c and at the back for F > c, as opposed to the 2d 
system which exhibits compression at the back for V < c and 
front for F > c (Fig. [^. 


riod 2iIF, of image dipoles located at z„ = Xn + it/™ 
and Zn = Xn + — yn), respectively (Fig. &)■ To cal¬ 

culate the velocity held w due to N swimmers and their 
image system, one has to evaluate inhnite sums of terms 
that decay as 1 jz^. An approximate numerical solution is 
given in [20] . Here, we derive an exact analytical solution 
in terms of hyperbolic functions. 
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where aj = t7e‘“j — (1 — p)V and a^aj indicates the 
strength and orientation of the dipole induced at Zj. 

We examine the emergent behavior of populations of N 
swimmers that initially hll the whole channel width and 
a segment £ of the channel length and are homogeneously 
placed and randomly orientated. Periodic boundary con¬ 
ditions are imposed at the two ends of the channel. We 
set p — 0.5. When the swimmers are subject to a suf- 
hciently strong background how, they quickly reorient 
and align in either the same (for z/ > 0) or opposite (for 
z/ < 0) direction to the background flow. Basically, the 
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background flow suppresses the orientation dynamics of 
the swimmers as noted in |20j . All aj reach, in finite time, 
a developed value ad = pU — (1 — p)V, where p = sgn(i/) 
and sgn is the signup function. It turns out that Cd, 
which we refer to as the effective polarization parame¬ 
ter, is a key parameter for understanding the emergent 
structure of the swimmers. 

The swimmers self-organize into a shock wave struc¬ 
ture where, for v > 0, the location of the compression 
wave depends on the strength of the background flow V, 
which in turn dictates the ‘speed of sound’ of the medium 
given by c = U + pV (see Supp. Mat. and Supp. Movie 
1). For a weak background flow V < c, the effective po¬ 
larization ad of the swimmers is dominated by their self- 
propelled motion so that ad > 0. The emergent structure 
is characterized by a compression wave at the back and 
a rarefaction or expansion wave at the front (Fig. [2 r). 
The average population speed is smaller than c, thus this 
shock wave can be termed as ‘subsonic’. If we tune the 
background flow properly such that the polarizations in¬ 
duced by the self-propelled motion and background flow 
cancel out such that ad = 0, the shock wave is suppressed 
and the swimmers evolve as a uniform cluster traveling 
at the speed V (Fig. [^). When the background flow is 
strong > c, it dominates the swimmers’ effective polar¬ 
ization ad < 0; the density shock forms at the front rather 
than at the back (Fig. and the average population 
speed is greater than c, thus the term ‘supersonic.’ Note 
that for ly < 0, the swimmers align in the opposite direc¬ 
tion to the background flow. The developed polarization 
parameter ad is always negative, and thus the swimmers 
develop a shock wave structure similar to Fig. i: that 
could travels upstream or downstream depending on the 
value of V. 

The density shock formed in Fig. is reminiscent to 
the experimental results of [3] for driven droplets in a 
Hele-Shaw channel and of the propagating density bands 
shown in supp. videos] for self-propelled colloids. 
The propagating density fronts in the latter system and 
their similarity to our results are surprising; however, 
the agreement with [3] is consistent with the intuition 
that active swimmers behave like passive particles when 
subject to strong background flow. 

We now consider the simple picture of three swimmers 
perfectly aligned parallel to the x-axis. One can easily 
see that for 1^ < c and ad > 0, HI push the swimmers 
forward such that the middle swimmer is always faster 
leading to compression at the front, whereas for V > c 
and ad < 0, HI hinder the swimming motion and the 
middle swimmer is always slower leading to compression 
at the back (Fig. [^. This argument can be generalized 
to N swimmers. Thus, in Id systems, the density shock 
forms at the front of the swimmers when the background 
flow is weak and at the back when the background flow 
is strong. This prediction, illustrated in Fig. agrees 
with the density shocks seen in the Id experiments of 
driven microfiuidic droplets [33], and is exactly opposite 
to the 2d experiments [3] and simulations presented in 


Fig. [2] The contradiction reveals the importance of 2d 
HI and the complexity of the mechanism responsible for 
the formation of density shock waves in populations of 
microswimmers in confined 2d channels. 

Density shocks emerge as a result of the interplay be¬ 
tween 2d hydrodynamics, sidewall confinement, and the 
initial distribution of swimmers into a segment i of the 
full channel. For W = oo and swimmers homogeneously 
filling the whole space, no shock develops due to the 
radial symmetry of the dipolar field at each swimmer. 
When sidewall confinement is introduced, it serves as a 
mechanism to break the radial symmetry of the dipo¬ 
lar field and creates a non-zero average hydrodynamic 
interaction acting on the swimmers. The initial uniform 
distribution of swimmers into a segment i of the full chan¬ 
nel length creates a spatial gradient in the HI in the x- 
direction, which contributes to triggering the instability 
of the initially uniform swimmer distribution. 

We numerically investigate the effect of 2d interactions 
on the shock wave formation by fixing = 1 and V = 1 
(i.e., ad = 0.5) and varying the channel width W and 
the initial local area fraction = ira^N/^Wi). To ob¬ 
tain the desired value of we set £ = 100 and vary 
N. For each set of parameters, we run multiple trials 
corresponding to different sets of initial conditions taken 
from a uniform probability distribution function (Monte- 
Carlo type simulations). We capture the emergence of 
the density shock wave using the shock order parameter 

S{t) = 1 -2'^\p{x,t) - pfit{x,t)\ /'^pfit{x,t) (3) 

X X 

where p(x) is the local density of the swimmers and 
Pfd.{x) is a triangular fit of the density profile p{x) 
(Fig. 1^). S' = 1 corresponds to a perfectly triangular 
shock and S = 0 to a rectangular homogenous profile. An 
example of the evolution of S is depicted in Fig. j^. The 
dependence of the mean value {S)t^oo averaged over all 
simulations on the parameter space (VF, d)^) is depicted 
in Fig. [^. It shows a continuous increase in (S)t_>oo as 
W and increase which then levels off for large W and 
This plot provides two important pieces of informa¬ 
tion. First, it serves as a strong supporting evidence that 
2d HI is important for the shock wave formation - for a 
given an ordered shock is only formed when W ex¬ 
ceeds a critical value. Second, it shows that the density 
shock formation is robust over a large plateau of W and 

To further elucidate the effect of the complex inter¬ 
play between the 2d HI, sidewall confinement, and back¬ 
ground flow on the density shock formation, we develop 
a reduced continuum model that properly accounts for 
all these elements. Our model is distinct from the traffic 
flow model used in and Burgers’ equation used in [3]. 
These equations oversimplify the effect of HI. They are 
only correct to leading order in the weak confinement 
limit with a W and are not capable of explaining 
the transition of shock formation from subsonic to super- 
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FIG. 4. (a) Shock order parameter is defined by comparing 
the actual density profile to a triangular density profile, shown 
here at t = 2000 for IF = 80 and = 0.3. (b) Evolution of 
the shock order parameter for the case in (a), (c) Mean value 
of the shock order parameter as a function of channel width 
W and swimmers’ area fraction 4 >a. 


sonic as the background flow velocity V increases. Here, 
we assume the swimmers’ density is homogenous in the 
y-direction and approximate the swimmers’ number den¬ 
sity as a Id function pc{x,t). The continuum number 
density pc is related to p via pc = pf{'iTa^). Following 
a standard homogenization procedure, we develop a Id 
continuum model (recall that p = sgn(:/)) 


9pc 

dt 


dx 


{Up + pV + pu) Pc 


= 0. 


(4) 


The 2d HI are captured by the term u{x, t) 


i{x, t) = 


1 


W-2a 


-Re 


Wc{z, z')pc{x', t)dydy'dx' 


(5) 


where x' is integrated from 0 to L and y, y' from —W/2+a 
to IF/2 — a, while the dipolar flow field Wc is given by 



We account for the steric repulsion between the swim¬ 
mers via the excluded area \x — x'\^ + \y — y'^ < 4a^, 
inside which the dipolar contribution is removed. The 


effect of the excluded area is also reflected in the up¬ 
per and lower limits of the y-integral in ([^, indicating 
that no swimmers can penetrate the sidewall boundaries. 
We solve the continuum equation 0 numerically with 
periodic boundary conditions. The numerical solutions, 
superimposed in red in Figure show excellent agree¬ 
ment between the discrete particle simulations and the 
continuum model, including successfully capturing the 
transition in the location of the shock formation as the 
background flow increases and the time scale of wave dis¬ 
persion (see Supp. Movie 1). 

This letter provides the first systematic study of the 
emergence of compression shock waves in populations of 
mircoswimmers driven in narrow flow channels. Our re¬ 
sults are consistent with experimental observations on 
driven droplets and self-propelled colloids, [siisiiiHiiin], 
but go beyond these observations to report an unob¬ 
served transition from ‘subsonic’ to ‘supersonic’ compres¬ 
sion shocks as the intensity of the background flow in¬ 
creases. The physical mechanisms underlying this tran¬ 
sition are elucidated via discrete simulations and a con¬ 
tinuum model that properly captures the effects of HI and 
geometric confinement. A rigorous analysis of the onset 
and propagation of instability in the continuum model is 
deferred to future work. 

This physical understanding of the emergent behavior 
and its dependence on the swimmers’ motility proper¬ 
ties and intensity of background flow can be exploited to 
devise novel mechanisms for controlling populations of 
microswimmers in externally-driven flow channels. Two 
species of microswimmers of different motility properties 
subject to the same background flow V get segregated in 
space in hnite time (see Supp. Movie 2). This is partic¬ 
ularly useful for biomedical applications such as sorting 
of cells in flow channels. Another compelling scenario is 
to control the throughput, i.e., the speed of a cell popu¬ 
lation through a flow channel, by careful manipulation of 
the background flow V. One could also control the pro¬ 
file of the population density distribution. For example, 
a Gaussian density profile can be obtained and main¬ 
tained from an initially uniform distribution by properly 
tailoring a time-dependent background flow V (see Supp. 
Movie 3). These results, albeit in the context of a simpli¬ 
fied model, have profound implications on developing a 
physics-based, quantitative framework for the design and 
control of particle-laden flows in microfluidic channels. 
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